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Mixed-mode oscillations and interspike interval statistics in 
the stochastic FitzHugh-Nagumo model 

Nils BergluncK and Damien Landon*' 



. Abstract 



We study the stochastic FitzHugh-Nagumo equations, modelling the dynamics of 
neuronal action potentials, in parameter regimes characterised by mixed-mode oscilla- 
tions. The interspike time interval is related to the random number of small-amplitude 
lO \ oscillations separating consecutive spikes. We prove that this number has an asymp- 

totically geometric distribution, whose parameter is related to the principal eigenvalue 
of a substochastic Markov chain. We provide rigorous bounds on this eigenvalue in 



the small-noise regime, and derive an approximation of its dependence on the system's 
parameters for a large range of noise intensities. This yields a precise description of 
the probability distribution of observed mixed-mode patterns and interspike intervals. 
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1 Introduction 

Deterministic conduction-based models for action-potential generation in neuron axons 
have been much studied for over half a century. In particular, the four-dimensional 
Hodgkin-Huxley equations [HH52J have been extremely successful in reproducing the ob- 
served behaviour. Of particular interest is the so-called excitable regime, when the neuron 
is at rest, but reacts sensitively and reliably to small external perturbations, by emitting a 
so-called spike. Much research efforts have been concerned with the effect of deterministic 
perturbations, though the inclusion of random perturbations in the form of Gaussian noise 
goes back at least to [GM64j . A detailed account of different models for stochastic per- 
turbations and their effect on single neurons can be found in [Tuc89] . Characterising the 
influence of noise on the spiking behaviour amounts to solving a stochastic first-exit prob- 
lem [Tuc75] . Such problems are relatively well understood in dimension one, in particular 
for the Ornstein-Uhlenbeck process |CR71| lTuc771 IRS80J . In higher dimensions, however, 
the situation is much more involved, and complicated patterns of spikes can appear. See 
for instance |TP01b| ITTP02I IRow07| for numerical studies of the effect of noise on the 
interspike interval distribution in the Hodgkin-Huxley equations. 
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FIGURE 1. Schematic phase diagram of the stochastic FitzHugh-Nagumo equations. The 
parameter a measures the noise intensity, S measures the distance to the singular Hopf 
bifurcation, and e is the timescale separation. The three main regimes are characterised 
be rare isolated spikes, clusters of spikes, and repeated spikes. 



Being four-dimensional, the Hodgkin-Huxley equations are notoriously difficult to 
study already in the deterministic case. For this reason, several simplified models have 
been introduced. In particular, the two-dimensional FitzHugh-Nagumo equations [Fit55^ 
Fit6lj lNAY62j . which generalise the Van der Pol equations, are able to reproduce one type 
of excitability, which is associated with a Hopf bifurcation (excitability of type II [IzhOOj). 

The effect of noise on the FitzHugh-Nagumo equations or similar excitable systems has 
been studied numerically [Lon931 IKP031 IKP061 ITGOS081 IBKLLClTj and using approxi- 
mations based on the Fokker-Planck equations [LSG99, SKllj . moment methods [T POlal 
TRW03], and the Kramers rate |LonOO| . Rigorous results on the oscillatory (as opposed to 
excitable) regime have been obtained using the theory of large deviations [MVEE05 , DT09] 
and by a detailed description of sample paths near so-called canard solutions [Sow08j. 

An interesting connection between excitability and mixed-mode oscillations (MMOs) 
was observed by Kosmidis and Pakdaman [KP03, KP06], and further analysed by Mura- 
tov and Vanden-Eijnden [MVE08] . MMOs are patterns of alternating large- and small- 
amplitude oscillations (SAOs), which occur in a variety of chemical and biological systems 
[DOP791 IHHM791 IPSS921 iDMS+OOj ■ In the de terministic case, at least three variables are 
necessary to reproduce such a behaviour (see |DGK + lT for a recent review of determin- 
istic mechanisms responsible for MMOs) . As observed in [KP031 IKP061 IMVE08| , in the 
presence of noise, already the two-dimensional FitzHugh-Nagumo equations can display 
MMOs. In fact, depending on the three parameters noise intensity a, timescale separation 
e and distance to the Hopf bifurcation 5, a large variety of behaviours can be observed, 
including sporadic single spikes, clusters of spikes, bursting relaxation oscillations and co- 
herence resonance. Figure [1] shows a simplified version of the phase diagram proposed 
in [MVE08| . 

In the present work, we build on ideas of |MVE08j to study in more detail the transition 
from rare individual spikes, through clusters of spikes and all the way to bursting relaxation 
oscillations. We begin by giving a precise mathematical definition of a random variable N 
counting the number of SAOs between successive spikes. It is related to a substochastic 
continuous-space Markov chain, keeping track of the amplitude of each SAO. We use this 
Markov process to prove that the distribution of N is asymptotically geometric, with a 
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parameter directly related to the principal eigenvalue of the Markov chain (Theorem 13, 2\i . 
A similar behaviour has been obtained for the length of bursting relaxation oscillations in 
a three-dimensional system [HM09J . In the weak-noise regime, we derive rigorous bounds 
on the principal eigenvalue and on the expected number of SAOs (Theorem 14. 2p . Finally, 
we derive an approximate expression for the distribution of N for all noise intensities up 
to the regime of repeated spiking (Proposition 15. 

The remainder of this paper is organised as follows. Section [2] contains the precise 
definition of the model. In Section O we define the random variable N and derive its 
general properties. Section 0] discusses the weak-noise regime, and Section [5] the transition 
from weak to strong noise. We present some numerical simulations in Section [6j and give 
concluding remarks in Section [7J A number of more technical computations are contained 
in the appendix. 
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We will consider random perturbations of the deterministic FitzHugh-Nagumo equations 



where a,b,c E K and e > is a small parameter. The smallness of e implies that x changes 
rapidly, unless the state (x,y) is close to the nullcline {y = x 3 — x}. Thus System (|2.ip is 
called a fast-slow system, x being the fast variable and y the slow one. 

We will assume that b ^ 0. Scaling time by a factor b and redefining the constants a, c 
and e, we can and will replace b by 1 in f)2 . 1 [) . If c ^ and c is not too large, the nullclines 
{y = x 3 — x} and {a = x + cy} intersect in a unique stationary point P. If c < 0, the 
nullclines intersect in 3 aligned points, and we let P be the point in the middle. It can be 
written P = (a, a 3 — a), where a satisfies the relation 
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given by 



ex = x — x 3 + y 
y = a — bx — cy 



(2.1) 



a + c(a 3 — a) 



= a . 



(2.2) 



The Jacobian matrix of the vector field at P is given by 




(2.3) 
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Figure 2. Some orbits of the deterministic FitzHugh-Nagumo equations (|2.1| for param- 
eter values e — 0.05, a — 0.58, 6 = 1 and c = 0. The black curve is the nullcline, and the 
red orbit is the separatrix. 

It has determinant (1 — c(l — 3a 2 )) je and trace 



Tr J = v * g - , where a* = W — - — . (2.4) 

Thus if |c| < \/-*Je, ^ admits a pair of conjugate imaginary eigenvalues when a = ±a». 
Furthermore, the eigenvalues' real parts are of order (a* — a)/e near a*. The system 
undergoes so-called singular Hopf bifurcations [BE86, B E921 lBra98| at a = ±o#. 

We are interested in the excitable regime, when a — a* is small and positive. In this 
situation, P is a stable stationary point, corresponding to a quiescent neuron. However, a 
small perturbation of the initial condition, e.g. a slight decrease of the ^-coordinate, causes 
the system to make a large excursion to the region of negative x, before returning to P 
(Figure [2]). This behaviour corresponds to a spike in the neuron's membrane potential, 
followed by a return to the quiescent state. One can check from the expression of the 
Jacobian matrix that P is a focus for a — a* of order \fe. Then return to rest involves 
small-amplitude oscillations (SAOs), of exponentially decaying amplitude. 

For later use, let us fix a particular orbit delimiting the spiking and quiescent regimes, 
called separatrix. An arbitrary but convenient choice for the separatrix is the negative- 
time orbit of the local maximum 

(-l/\/3,2/(3\/3)) of the nullcline (Figure [2]). The main 
results will not depend on the detailed choice of the separatrix. 

In this work we consider random perturbations of the deterministic system (|2.ip by 
Gaussian white noise. They are described by the system of ltd stochastic differential 
equations (SDEs) 

dx t = -(x t -x 3 t + y t ) dt + ^= dW t {1) 

£ (2.5) 

(2) 

dy t = (o - x t - cy t ) dt + a 2 dWj: , 

where and are independent, standard Wiener processes, and o"i,o"2 > 0. The 

parameter a will be our bifurcation parameter, while c is assumed to be fixed, and small 
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Figure 3. Examples of time series of the stochastic FitzHugh-Nagumo equations ()2.5[) . 
The plots show the functions t H> £t, where the variable £ is defined in Section^ Parameter 
values are e = 0.01 and 5 = 3- 10~ 3 for the top row, 5 — 5 ■ 10~ 3 for the bottom row. The 
noise intensities are given by o\ = a 2 = 1.46 • 10~ 4 , 1.82 • 10~ 4 , 2.73 ■ 10~ 4 and 3.65 • 10~ 4 . 



enough for the system to operate in the excitable regime. The scaling in 1 / yfe of the noise 
intensity in the first equation is chosen because the variance of the noise term then grows 
like crft/e, so that a\ measures the ratio of diffusion and drift for the x-variable, while <t| 
plays the same role for the y- variable. 

Figure [3] shows a selection of time series for the stochastic FitzHugh-Nagumo equa- 
tions (|2.5p . For the chosen parameter values, one can clearly see large-amplitude spikes, 
separated by a random number of SAOs. Note the rather large variability of the SAOs' 
amplitude. 



3 The distribution of small-amplitude oscillations 

In this section we define and analyse general properties of an integer-valued random vari- 
able N, counting the number of small-amplitude oscillations the stochastic system per- 
forms between two consecutive spikes. The definition is going to be topological, making 
our results robust to changes in details of the definition. We start by fixing a bounded 
set T> C M 2 , with smooth boundary &D, containing the stationary point P and a piece 
of the separatrix (Figure 0J). Any excursion of the sample path (xt,yt)t outside V will be 
considered 

To define precisely, we let B be a small ball of radius p > centred in P. Then we 
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Figure 4. Definition of the number N of SAOs. The sample path (blue) enters the region 
T>, and intersects twice the line J- before leaving T>, making another spike. Thus N = 2 
in this example. The separatrix is represented in red. 

draw a smooth curve T from B to the boundary dT>, which we parametrise by a variable 
r G [0, 1] proportional to arclength (the results will be independent, however, of the choice 
of T and of r). We extend the parametrisation of J 7 to a polar-like parametrisation of 
all T>\B, i.e. we choose a diffeomorphism T : [0, 1] x § 1 — > V, (r,ip) i— > (x,y), where 
T~ l (F) = {ip = 0}, T~ 1 (&D) = {r = 0} and T~ x (dB) = {r = 1}. We also arrange that 
if > near P for the deterministic flow. 

Consider the process (rt, <pt)t (where the angle <p has been lifted from S 1 to M). Given 
an initial condition (ro, 0) G T _1 (J r ) and an integer M ^ 1, we define the stopping time 

r = inf {t > 0: <p t € {2vr, -2Mvr} or r t G {0, 1}} . (3.1) 

There are four cases to consider: 

• The case r T = corresponds to the sample path (xt,yt) leaving V, and thus to a 
spike. This happens with strictly positive probability, by ellipticity of the diffusion 
process (12. 5h . In this situation, we set by convention N = 1. 

• In the case ip T = 2ir and r T G (0, 1), the sample path has returned to J- after performing 
a complete revolution around P, staying all the while in T> \B. This corresponds to 
an SAO, and thus N ^ 2. 

• The case r T = 1 corresponds to the sample path entering B, which we consider as the 
neuron reaching the quiescent state. In that case we simply wait until the state leaves 
B again and either hits J- or leaves T>. 

• The case <p T = —2Mir and r T G (0, 1) represents the (unlikely) event that the sample 
path winds M time around P in the wrong direction. We introduce this case for 
technical reasons only, as we will need r to be the first-exit time of a bounded set. For 
simplicity, we also consider this situation as one SAO. 

As long as r T G (0, 1), we repeat the above procedure, incrementing N at each iteration. 
This yields a sequence (Ro, R±, . . . , Rn-i) of random variables, describing the position of 
the successive intersections of the path with P, separated by rotations around P, and up 
to the first exit from T>. 
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Remark 3.1. The above definition of N is the simplest one to analyse mathematically. 
There are several possible alternatives. One can, for instance, introduce a quiescent state 
(x, y) G £>, and define N as the number of SAOs until the path either leaves T> or enters B. 
This would allow to keep track of the number of SAOs between successive spikes and/or 
quiescent phases. Another possibility would be to count rotations in both the positive and 
negative directions. For simplicity, we stick here to the above simplest definition of N, 
but we plan to make a more refined study in a future work. 

The sequence (Rn)n forms a substochastic Markov chain on E — (0, 1), with kernel 

K(R, A) = ¥{(p T =G {2vr, -2Mvr}, r T E A \ (p = 0, r = R} , 

R G E, A C E a Borel set . (3.2) 

The Markov chain is substochastic because K{R,E) < 1, due to the positive probability 
of sample paths leaving T>. We can make it stochastic in the usual way by adding a 
cemetery state A to E (the spike), and setting K(R, A) = 1 - K(R,E), K(A, A) = 1 
(see [Ore714 lNum84| for the general theory of such processes). 
The number of SAOs is given by 

N = inf {n > 1 : R n = A} G N U {oo} (3.3) 

(we set inf = oo). A suitable extension of the well-known Perron-Frobenius theorem (see 
[Jenl21 IKR501 IBir57] ) shows that K admits a maximal eigenvalue Ao , which is real and 
simple. It is called the principal eigenvalue of K. If there exists a probability measure 
ttq such that ttqK = Xottq, it is called the quasi- stationary distribution (QSD) of the 
kernel K |SV,T66j . 

Our first main result gives qualitative properties of the distribution of TV valid in all 
parameter regimes with nonzero noise. 

Theorem 3.2 (General properties of N). Assume that 01,0-2 > 0. Then for any initial 
distribution fiQ of Rq on the curve T , 

• the kernel K admits a quasi- stationary distribution ttq; 

• the associated principal eigenvalue Ao = Ao(e, a, c, 0*1, 02) is strictly smaller than 1; 

• the random variable N is almost surely finite; 

• the distribution of N is "asymptotically geometric", that is, 

lim P w, {iV = n + l\N > n\ = 1 - A ; (3.4) 

• E^ {r^} < 00 for r < 1/Xq and thus all moments E^ |iV fe } of N are finite. 

Proof: Let us denote by K(x,dy) the kernel defined in (13.2p . We consider if as a 
bounded linear operator on L°°(E), acting on bounded measurable functions by 

f(x) 1 ^ (Kf)(x) = [ K(x,dy)f(y)=E x {f(R 1 )} , (3.5) 

and as a bounded linear operator on L l {E), acting on finite measures fi by 

fi(A) QiK)(A) = [ n(dx)K(x, A) = ¥^{R 1 G A} . (3.6) 
Je 
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To prove existence of a QSD ttq, we first have to establish a uniform positivity condition 
on the kernel. Note that in K(x,dy), y represents the first-exit location from the domain 
Q = (0, 1) x (-2M-K, 2tt), for an initial condition (x, 0), in case the exit occurs through one 
of the lines <p = —2Mir or (p = 2ir. In harmonic analysis, K(x, dy) is called the harmonic 
measure for the generator of the diffusion in Q based at (x,0). In the case of Brownian 
motion, it has been proved in |Dah77] that sets of positive Hausdorff measure have positive 
harmonic measure. This result has been substantially extended in [BAKS84 , where the 
authors prove that for a general class of hypoelliptic diffusions, the harmonic measure 
admits a smooth density k(x, y) with respect to Lebesgue measure dy. Our diffusion 
process being uniformly elliptic for a\ y ai > 0, it enters into the class of processes studied 
in that work. Specifically, |BAKS84~1 Corollary 2.11] shows that k(x,y) is smooth, and its 
derivatives are bounded by a function of the distance from x to y. This distance being 
uniformly bounded below by a positive constant in our setting, there exists a constant 
L G R + such that 

sup k(x, y) 

; " /: - . /. VxeE. (3.7) 



inf k(x, y) 



We set 



s(x) = inf k(x,y) . (3-8) 

yeE 

Then it follows that 

s(x) ^ k(x,y) ^ Ls(x) Vx,y G E . (3-9) 
Thus the kernel K fulfils the uniform positivity condition 

s(x)v(A) ^ K(x,A) ^ Ls{x)v(A) \/x G E ,VA C E (3.10) 

for v given by the Lebesgue measure. It follows by |Bir57l Theorem 3] that K admits 
unique positive left and right unit eigenvectors, and that the corresponding eigenvalue Ao 
is real and positive. In other words, there is a measure ttq and a positive function ho such 
that ttqK = Aqvto and Kho = Aq/io- We normalise the eigenvectors in such a way that 



ir (E) = / 7r (d2;) = 1 , Tr h = / ir (dx)h (x) = 1 . (3.11) 
JE JE 

Thus 7To is indeed the quasistationary distribution of the Markov chain. Notice that 

A = XoME) = [ ir (dx)K(x,E) < ir (E) = 1 , (3.12) 

JE 

with equality holding if and only if K(x, E) = 1 for 7To-almost all x G E. In our case, 
K(x,E) < 1 since E has strictly smaller Lebesgue measure than dQ, and the density of 
the harmonic measure is bounded below. This proves that Ao < 1. 

We denote by K n the n-step transition kernel, defined inductively by K 1 = K and 

K n+1 (x,A)= [ K n (x,dy)K(y,A) . (3.13) 

JE 

Lemma 3 in [Bir57j shows that for any bounded measurable function / : E — > R , there 
exists a finite constant M(f) such that the spectral-gap estimate 

\(K n f)(x) - KM)h (x)\ < M(f)(\ p) n h (x) (3.14) 
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holds for some p < 1 (note that this confirms that Ao is indeed the leading eigenvalue of 
K). In order to prove that A is almost surely finite, we first note that 

P^°{A >n}= F^ {R n £E} = [ no(dx)K n (x,E) . (3.15) 

Je 

Applying (|3.14p with / = 1, the function identically equal to 1, we obtain 

ASM*) - M(l)(X oP ) n h (x) ^ K n (x, E) ^ \%h (x) + M(l)(\ oP ) n h {x) • (3.16) 

Integrating against no, we get 

A£(l - M{l)p n ) < P«'{A > n) < A£(l + M{l)p n ) . (3.17) 

Since Ao < 1, it follows that lim^—^oo F flo {N > n} = 0, i.e., A" is almost surely finite. 
In order to prove that A" is asymptotically geometric, we have to control 

F^{N = n + l}=[ [ p (dx)K n (x,dy)[l- K(y,E)] . (3.18) 
Je Je 

Applying (I3.14D with f(y) = 1 — K(y,E), and using the fact that 

nof = 1 - ( Mdy)K(y, E) = 1 - A (3.19) 
Je 

yields 

Ag(l - Ao - M(/)p n ) < P«{A = n + 1} < Ag(l - A + M(/)p") . (3.20) 

Hence (|3.4p follows upon dividing (|3.20p by (|3.17p and taking the limit n — >• oo. 

Finally, the moment generating function E w {r^} can be represented as follows: 



E^°{r N } = Y j r n ^°{N = n} = 1 + (r - 1) 



n-1 



m=0 



= 1 + (r- 1) ^r m P w {A > to} , (3.21) 

which converges for |rAo| < 1 as a consequence of (I3.17p . □ 

Note that in the particular case where the initial distribution po is equal to the QSD 
7To, the random variable R n has the law p n = Ag7ro, and A follows an exponential law of 
parameter 1 — Ao : 

¥ 7ro {N = n} = AS -1 (l-Ao) and E W °{N} = — — . (3.22) 

1 — Ao 

In general, however, the initial distribution po after a spike will be far from the QSD ttq, 
and thus the distribution of N will only be asymptotically geometric. 

Theorem 13.21 allows to quantify the clusters of spikes observed in [MVE08] . To this 
end, we have to agree on a definition of clusters of spikes. One may decide that a cluster 
is a sequence of successive spikes between which there is no complete SAO, i.e. A = 1 
between consecutive spikes. If the time resolution is not very good, however, one may 
also fix a threshold SAO number no ^ 1, and consider as a cluster a succession of spikes 
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separated by at most no SAOs. Let ^ be the arrival distribution on J- of sample paths 
after the nth spike. Then the probability to observe a cluster of length k is given by 

(o) , (l) , , O-i) , (fc) , 

P^o {iV^no}P^ {N < n } ...P M o {N ^ n }P M ° {N > n ] . (3.23) 

(n) 

In general, the consecutive spikes will not be independent, and thus the distributions Hq 
will be different. For small noise, however, after a spike sample paths strongly concentrate 
near the stable branch of the nullcline (see the discussion in |BG09l Section 3.5.2]), and 

(n) 

thus we expect all fi to be very close to some constant distribution fiQ. This implies 
that the lengths of clusters of spikes also follow an approximately geometric distribution : 

Pjcluster of length k] ~ p k (l - p) where p = P M {iV ^ n } . (3.24) 



4 The weak-noise regime 

In order to obtain more quantitative results, we start by transforming the FitzHugh- 
Nagumo equations to a more suitable form. The important part of dynamics occurs near 
the singular Hopf bifurcation point. We carry out the transformation in four steps, the 
first two of which have already been used in [BE861 |BE92 : 

1. An affine transformation x = a* + u, y = a* — a* + v translates the origin to the 
bifurcation point, and yields, in the deterministic case (|2.1|) . the system 

eu = v + ceu — 3a* -u — u , 

(4.1) 

v = 5 — u — cv , 

where S = a — a* — c(a* — a*) is small and positive. Note that (|2.2p implies that 5 is 
of order a — a* near the bifurcation point, and thus measures the distance to the Hopf 
bifurcation. In particular, by (|2.4j) the eigenvalues of the Jacobian matrix J have real 
parts of order —5/e. 

2. The scaling of space and time given by u = y/e£, v = er\ and t = yfet' yields 

£ = ??-3a*£ 2 + ^(c£-£ 3 ) , 

• 6 f r (42) 
n = ~r - £ - V £ cr] , 

v e 

where dots now indicate derivation with respect to t' . On this scale, the nullcline £ = 
is close to the parabola r] = 3a* £ 2 . 

3. The nonlinear transformation rj = 3a* £ 2 + z — l/(6a*) has the effect of straightening 
out the nullcline, and transforms (14.21) into 



(4.3) 



z = — p - 6a*£z + \fe{ 6a*£ 4 + c(— ^ 9a*£ 2 

\/e \ V6a* 

4. Finally, we apply the scaling £ h-> — £/3a*, z \-t z/3a*, which yields 



9a 2 

, + 2e, + ^(^ 4 + c (i-3e 2 



(4.4) 
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Figure 5. (a) Level curves of the first integral Q = 2ze~ 2z ~ 2 ^ +1 . (b) Some orbits 
of the deterministic equations (|4.4p in (£, z)-coordinates, for parameter values e = 0.01, 
|U = 0.05 and c = (i.e. a* = l/y/3). 



where the distance to the Hopf bifurcation is now measured by the parameter 

3a *^ fA r\ 

M = ■ (4.5) 

Let us first consider some special cases of the system (|4.4p . 

• If e = // = 0, we obtain 

5 2 (4.6) 
z = 2£z . 

This system admits a first integral 

Q = 2ze- 2z - 2 f +1 , (4.7) 

which is equivalent to the first integral found in [B E92] (we have chosen the normal- 
isation in such a way that Q £ [0,1] for z ^ 0). Figured shows level curves of Q. 
The level curve Q = corresponds to the horizontal z = 0. Values of Q 6 (0, 1) yield 
periodic orbits, contained in the upper half plane and encircling the stationary point 
P = (0,1/2) which corresponds to Q = 1. Negative values of Q yield unbounded 
orbits. Hence the horizontal z = acts as the separatrix in these coordinates. 

• If \i > and e = 0, the stationary point P moves to (—(i, 1/2) and becomes a focus. 
The separatrix is deformed and now lies in the negative half plane. It delimits the 
basin of attraction of P. 

• If \x > and < £ < 1, the dynamics does not change locally (Figure [5b). The global 
dynamics, however, is topologically the same as in original variables. Therefore, orbits 
below the separatrix are no longer unbounded, but get "reinjected" from the left after 
making a large excursion in the plane, which corresponds to a spike. Orbits above the 
separatrix still converge to P, in an oscillatory way. 

Carrying out the same transformations for the stochastic system (j2.5j) yields the fol- 
lowing result (we omit the proof, which is a straightforward application of Ito's formula). 
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FIGURE 6. Sample paths of the stochastic equations (|4.8[) (black), superimposed on the 
level curves of Q (blue). The red line is the separatrix for the case e = p, = 0. Parameter 
values are c = 0, e — 0.01 and (a) p — 0.01, o\ = a% = 0.01, (b) p = 0.01, ax = a% = 0.03, 
(c) p = 0.05, ax = a 2 = 0.01, (d) p = -0.05, a x = a 2 = 0.013. 



Proposition 4.1. In the new variables (£, z), and on the new timescale t/y/e, the stochas- 
tic FitzHugh-Nagumo equations (12, 5j) take the form 



dz t = 
where 



dt + cri dW, 



(i) 



/2 + 2^ + ^(^ 4 + c(i-3e, 2 



-3/4, 



2f 



dt - 2cri& dW t (1) + ct 2 dW t (2) 



(4. 



5-2 = 3a*£~ 3//4 er 2 , 

3a* (5 — 3a*<7 2 /£) 



(4.9) 



A» = M 



-2 



a 



Note the Ito-to-Stratonovich correction term — a\ in /2, which implies that this param- 
eter can be positive or negative, depending on the value of 5 and the noise intensity. 

Figure [6] shows sample paths of the SDE (|4.8p . superimposed on the level curves of 
the first integral Q. For sufficiently weak noise, the sample paths stay close to the level 
curves. Whether the system performs a spike or not depends strongly on the dynamics 
near the separatrix, which is close to z = when ft and e are small. 
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Figure 7. Definition of the curves F_ and F + . To approximate the kernel of the Markov 
chain, we take an initial condition zq on F_ , and determine first bounds on the first- hitting 
point Z\ of F + , and then on the first- hitting point z' of F_. 



To understand better the dynamics close to z = 0, consider the system 

1 



d6° = - dt , 

' 2 (4.10) 

OJTJ/W I ^^w( 2 ) 



dz t ° = fldt- 2a x § dW t { ' + o 2 dW t . 

obtained by neglecting terms of order z and yfe in (|4.8p . The solution of (|4.10p is given 
by 

rt (4.11) 
z° t =z% + jit- 2a! / g dW^ + a2 W t {2) . 
Jo 

There is a competition between two terms: the term jit, which pushes sample paths 
upwards to the region of SAOs, and the noise terms, whose variance grows like (af + cr|)i. 
We can thus expect that if \x > and af + a^ <C fl 2 , then the upwards drift dominates the 
standard deviation of the noise terms. The system will be steered to the upper half plane 
with high probability, and make many SAOs before ultimately escaping and performing a 
spike. Going back to original parameters, the condition a\ + a\ <C translates into 

al + al^ie 1 /^) 2 . (4.12) 

If by contrast a\ + a\ is of the same order as fi 2 or larger, or if jl < 0, then the probability 
of spiking will be much larger. 

We now make these ideas rigorous, by proving that spikes are indeed rare in the regime 
af + d\ <C jl 2 ■ In order to define the Markov chain, we introduce the broken line 

F_ = |^ = -L and z < ^} U j-L ^ £ < and z = -} , (4.13) 

see Figure [71 We choose an initial condition (— L,zq) on F_, with zq close to zero. Our 
aim is to estimate the probability that the sample path starting in (— L, zq) returns to F_, 
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after having made one revolution around P, at a point (— L, Z\) with Z\ > zq. This will 
be likely for weak noise. As a consequence, the set of points in i*L lying above (-L, zq) 
will be mapped into itself by the Markov chain with high probability, and this fact can be 
used to estimate the principal eigenvalue Ao and the expected number of SAOs. 

A difficulty of this approach is that the (£, z)-coordinates are only useful for small z 
and bounded £. For the remaining dynamics, it is in fact much simpler to work with 
the first integral Q (this idea has already been used in [MVE08I), and some conjugated 
angular variable (f>. It would be nice if the (Q, ^-coordinates could be used everywhere, 
but unfortunately it turns out they are not appropriate when z is close to (this is related 
to the fact that the period of oscillations diverges near the separatrix). We are thus forced 
to work with both pairs of variables, depending on where we are in the phase plane. So 
we introduce a second broken line 

F + = |^ = L and z s: ^} U jo s$ £ ^ L and z = ^} , (4.14) 

and use (£, z)-coordinates to describe the dynamics below F_ and F + , and ((^^-coor- 
dinates to describe the dynamics above these lines. 

An important point is the choice of the parameter L defining F±. On one hand, it 
cannot be too large, because we want to treat the terms of order \fe in (|4.8p as small 
perturbations. On the other hand, the equations in (Q, </>)-coordinates contain error terms 
which become large in the region of small z and bounded £. We will work with z bounded 
below by a positive constant cq times /i 1-7 for some appropriate 7 S (0, 1). Thus taking 

L 2 = |(-log(c^)) , (4.15) 

ensures that e -2 ^ 2 = (c_/i) 7 , and thus that the first integral Q has at least order fi when 
|£| ^ L and z ^ co/i 1-7 , where we will use (Q, ^-coordinates. 
We can now state the main result of this section. 

Theorem 4.2 (Weak-noise regime). Assume that e and 5/yfe are sufficiently small. If 
c / 0, assume further that 5 ^ \c\e^ for some f3 < 1. Then there exists a constant k > 
such that for a\ + a\ ^ (e l ^5) 2 / \og{y/e/5), the principal eigenvalue Ao satisfies 

1-Ao ^exp{- K (£ 2 1/4 ^ | . (4.16) 

Furthermore, for any initial distribution (j, of incoming sample paths, the expected number 
of SAOs satisfies 

E»°{N} > C(^ )expL (£ 2 1/ ^ ) 2 2 | . (4.17) 

Here C(/io) is the probability that the incoming path hits F_ above the separatrix. 

Proof: Let us first show that the problem can be reduced to proving the existence of a 
subset A C E with positive Lebesgue measure that the Markov chain is unlikely to leave. 
This set will then be chosen as the set of points in F- for which z is larger than some zq 
of order /i 1-7 . Let 

e A = sup[l - K(x,A)} (4.18) 
be the maximal probability to leave A when starting in A. Let us show that 

Xo^l-e A . (4.19) 
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Indeed, the relation Aovr = ir K yields 

Ao^o(^) = / n (dx)K(x,A) + / TT (dx)K(x,A) 

J A JE\A 

> 7T (A)(1 -e A )+ [ ir (dx)s(x) v{A) . (4.20) 
Je\a 

Either ttq(A) = 1, and the result follows immediately. Or ttq{A) < 1, and thus ttq(E\A) > 
0, so that the second term on the right-hand side is strictly positive. It follows that 
Ao7T"o(A) > 0, and we obtain (|4.19|) upon dividing by tto(A). 
Next, let us prove that 

E«> { N } ^ Mdl . (4.21) 

£A 

For x £ A, let 9{x) = E x {N} = ^ n>0 Rn ( x i E )- Then 9 ( x ) = lim n^oo 9 n (x) where 

n 

6 n {x) = Y,K m {x,E) . (4.22) 

m=0 

We have 

n+ i(a;) = 1 + (K9 n ){x) > 1+ / tf(a;,dy)0 n (y) . (4.23) 
Now let m n = inf^g^ 6 n (x). Then mo = 1 and 

m n+ i > 1 + (1 - e A )m n . (4.24) 

By induction on n we get 

m„ > J- - (1 ~ £aT+1 , (4.25) 
£A £A 

so that E 1 ' {iV} = 9{x) ^ for all x £ A, and (|4.21 j) follows upon integrating against 
Ho over A. 

It thus remains to construct a set A C E such that £a is exponentially small in 
jl 2 /(af + erf). The detailed computations being rather involved, we give them in the 
appendix, and only summarise the main steps here. 

1. In the first step, we take an initial condition (—L,zq) on i*L, with zq of order /i 1-7 
(Figure [7|). It is easy to show that the deterministic solution starting in (-L, zq) hits 
F + for the first time at a point Zj. ^ co/i 1-7 where cq > 0. Consider now the stochastic 
sample path starting in (— L, zq). Proposition IA.4I in Appendix lAl shows that there 
are constants C, k% > such that the sample path hits F + for the first time at a point 
(L, z\) satisfying 

¥{ Zl < COM 1-7 - A} < ^ e -^ 2 /° 2 . (4.26) 

This is done by first approximating (|4.8p by a linear system, and then showing that 
the effect of nonlinear terms is small. 

2. In the second step, we show that a sample path starting in (L, z{) G F + returns with 
high probability to _F_ at a point (— L, z' ) with z' ^ z\. Using an averaging procedure 
for the variables (Q, 4>), we show that Q varies little between F + and F_. Corollary IB. 51 
in Appendix IB1 shows that there is a Ki > such that 



{z' < Z1 \zx> cq/U 1-7 - ft} < 2e- K2 ^ 2/ ^ . (4.27) 
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In the above results, we assume that either c = 0, or c 7^ and ^fe for some 6 > 

(which follows from the assumption 5 > e^). The reason is that if c = 0, we can draw on 
the fact that the error terms of order yje in (|4.8p are positive, while if c 7^ we only know 
their order. Choosing A as the set of points in F_ for which z co/i 1 ^ 7 — fl, we obtain 
that £a is bounded by the sum of f|4.26j) and (|4.27p . and the results follow by returning 
to original parameters. □ 

Relation (|4.17|) shows that the average number of SAOs between two consecutive spikes 
is exponentially large in this regime. Note that each SAO requires a rescaled time of order 
1 (see Section lETTj) . and thus a time of order ^fe in original units. It follows that the average 
interspike interval length is obtained by multiplying (|4.17p by a constant times yfe. 

Relation (|3.4p shows that the distribution of N is asymptotically geometric with pa- 
rameter given by (|4,16p . Hence the interspike interval distribution will be close to an 
exponential one, but with a periodic modulation due to the SAOs. 

5 The transition from weak to strong noise 

We now give an approximate description of how the dynamics changes with increasing noise 
intensity. Assume that we start (I4.8D with an initial condition (£,o,zo) where £0 = —L for 
some L > and zq is small. As long as zt remains small, we may approximate £t in the 
mean by £0 + t/2, and thus z t will be close to the solution of 

dzj = (fi + tzj) dt - ait dW t W + (T 2 dW t {2) . (5.1) 



This linear equation can be solved explicitly. In particular, at time T = 4L, is close to 
L and we have the following result. 



Proposition 5.1. Let 2L 2 = 7|log(c_/i)| for some 7,c_ > 0. Then for any H, 



c {4 < -H) = $ 



-7T ' — 

a 



7-1 



(5.2) 



where a 2 = a\ + of and 3>(a?) = u ^ 2 du/y/2~K is the distribution function of the 

standard normal law. 



Proof: Solving (15. II) by variation of the constant yields 

f T e - 2 /2 ds - &1 [ T s e~ s2 / 2 dW^ + a 2 F , 

'to Jt() Jt() 



z\ = z + e T2/2 



(5.3) 



Note that by the choice of L, we have e T2//2 = e 2Z/2 = (c_/i) 7 . The random variable z\ 
is Gaussian, with expectation 

E{4} = Zo + fle 2L2 e" s2 / 2 ds (5.4) 

J-2L 

and variance 

Var(z^) = a\ e AL ' / s 2 e" 5 ' ds + ~a\ e iL * / e" 3 ' ds . (5.5) 

J-2L J-2L 
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Using this in the relation 



P{4 ^ --H"} = / ^=^= dz = $ 1 TS 

•/2vrVar(4) \ wVar(4) 



(5.6) 



yields the result. □ 



Choosing 7 large enough, the right-hand side of (|5.2p is approximately constant for a 
large range of values of zq and H. The probability that the system performs no complete 
SAO before spiking again should thus behave as 

F^{N = 1} * $ (-^1) = $ (" M ^|==^ M ) • (5-7) 

Since 1 — Ao is equal to the probability of leaving D before completing the first SAO, 
when starting in the QSD ttq, we expect that 1 — Ao has a similar behaviour, provided 
7i"o is concentrated near z = 0. We can identify three regimes, depending on the value of 
fr/v ■ 

1. Weak noise : jl S> <r, which in original variables translates into \J o\ + a 2 <C £ 
This is the weak-noise regime already studied in the previous section, in which Ao is 
exponentially close to 1, and thus spikes are separated by long sequences of SAOs. 

2. Strong noise : jl <C — cr, which implies fi <^ a 2 , and in original variables translates 
into yj o\ + a\ > e 3 / 4 . Then P{A > 1} is exponentially small, of order e~(. u l +a i)/ £3/2 . 
Thus with high probability, there will be no complete SAO between consecutive spikes, 
i.e., the neuron is spiking repeatedly. 

3. Intermediate noise : = 0(a), which translates into e l ^5 ^ V"i + a \ ^ £ 3 ^ 4 - 
Then the mean number of SAOs is of order 1. In particular, when <j\ = Ved, ji = 
and thus F{N = 1} is close to 1/2. 

An interesting point is that the transition from weak to strong noise is gradual, being 
characterised by a smooth change of the distribution of as a function of the parameters. 
There is no clear-cut transition at the parameter value o\ = \fs5 obtained in |MV E08j 
(cf. Figure [1]), the only particularity of this parameter value being that ¥{N = 1} is close 
to 1/2. In other words, the system decides between spiking and performing an additional 
SAO according to the result of a fair coin flip. The definition of a boundary between 
the intermediate and strong-noise regimes mainly depends on how well the SAOs can be 
resolved in time. A very good time resolution would put the boundary at noise intensities 
of order e 3 / 4 , while a lower time resolution would move in closer to ^fe5. 



6 Numerical simulations 

Figure [8] shows numerically simulated distributions of the SAO number. The geometric 
decay is clearly visible. In addition, for decreasing values of ju/(T, there is an increasing 
bias towards the first peak N = 1, which with our convention corresponds to the sys- 
tem performing no complete SAO between consecutive spikes. Of course this does not 
contradict the asymptotic result (13. 4p . but it shows that transient effects are important. 

Due to the finite sample size, the number of events in the tails of the histograms 
is too small to allow for a chi-squared adequacy test. We can, however, estimate the 
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Figure 8. Histograms of numerically simulated distributions of the SAO number N, 
obtained from time series containing 1000 spikes each. The superimposed curves show 
geometric distributions with parameter Ao, where Ao has been estimated from the expec- 
tation of r N , as explained in the text. Parameter values are e — 10~ 4 and a = 0.1 in all 
cases, and (a) p = 0.12, (b) fi = 0.05, (c) p = 0.01, and (d) P = -0.09 (cf. gSJ for their 
definition) . 



principal eigenvalue Ao, by using the fact that the moment generating function E^ {r^} 
has a simple pole at r = 1/Ao (see (|3.17|) and (|3.21|) ). Figure O shows examples of the 
dependence of the empirical expectation of r N on r. By detecting when its derivative 
exceeds a given threshold, one obtains an estimate of 1/Aq- Geometric distributions with 




Figure 9. Empirical expectation of r N as a function of r, for a = 0.1, e = 10 -4 and 
p = 0.05 (a), p — —0.03 (b) and p — —0.06 (c). The respective values of —p/a are thus 
—0.5 (a), 0.3 (b) and 0.6 (c). The location of the pole allows to estimate 1/Aq. 
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FIGURE 10. Plots of the function $(— i^'^ft/a) as a function of —ft/a (full line), and 
numerical estimates of ¥{N = 1} (stars), l/E{iV} (crosses) and 1 — A (circles). 

parameter Ao have been superimposed on two histograms in Figure El 

Figure [T0l shows, as a function of — x = —j2/a, the curve x h-> $(— vr 1 / 4 ^), as well as 
the inverse of the empirical expectation of N, the probability that N = 1, and 1 — Ao 
where the principal eigenvalue Ao has been estimated from the generating function. The 
data points for jj, > have been obtained from histograms containing 1000 spikes, while 
those for j2 < have been obtained from histograms containing 500 spikes separated by 
N > 1 SAOs (the number of spiking events with N = 1 being much larger). Theorem 14.21 
applies to the far left of the figure, when /i>(7. 

As predicted by (I5.7p . F{N = 1} is indeed close to the theoretical value <£(— Tr 1 / 4 /}/^). 
Recall from (|iT22]) that 1/E« {N}, V°{N = 1} and 1 - A would be equal if the initial 
distribution /j,q after a spike were equal to the QSD rro. The simulations show that l/E{iV} 
and 1 — Ao are systematically smaller than ¥{N = 1}. The difference between ¥{N = 1} 
and 1 — Ao is a measure of how far away fiQ is from the QSD ttq. The difference between 
l/E{iV} and 1 — Ao also depends on the spectral gap between Ao and the remaining 
spectrum of the Markov kernel. Note that l/E{A r } and 1 — Ao seem to follow a similar 
curve as ¥{N = 1}, but with a shifted value of fi. We do not have any explanation for 
this at the moment. 

7 Conclusion and outlook 

We have shown that in the excitable regime, and when the stationary point P is a fo- 
cus, the interspike interval statistics of the stochastic FitzHugh-Nagumo equations can 
be characterised in terms of the random number of SAOs N . The distribution of N is 
asymptotically geometric, with parameter 1 — Ao, where Ao is the principal eigenvalue 
of a substochastic Markov chain, describing a random Poincare map. This result is in 
fact fairly general, as it does not depend at all on the details of the system. It only re- 
quires the deterministic system to admit an invariant region where the dynamics involves 
(damped) oscillations, so that a Poincare section can be defined in a meaningful way. Thus 
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Theorem 13.21 will hold true for a large class of such systems. 

To be useful for applications, this qualitative result has to be complemented by quan- 
titative estimates of the relevant parameters. Theorem 14.21 provides such estimates for Ao 
and the expected number of SAOs in the weak- noise regime a\ + a\ <C {e l / A 5) 2 . We have 
obtained one-sided estimates on these quantities, which follow from the construction of an 
almost invariant region A for the Markov chain. It is possible to obtain two-sided estimates 
by deriving more precise properties for the Markov chain, in particular a lower bound on 
the probability of leaving the complement of A. We expect the exponent (e 1 / 4 ^) 2 / (c 2 +cr|) 
to be sharp in the case 5 ^> e, since this corresponds to the drift (I in the expression (|4.8p 
for z dominating the error terms of order y/e due to higher-order nonlinear terms. For 
smaller 5, however, there is a competition between the two terms, the effect of which is 
not clear and has to be investigated in more detail. The same problem prevents us from 
deriving any bounds for 5 ^ e when the parameter c defining the FitzHugh-Nagumo equa- 
tions is different from zero. It may be possible to achieve a better control on the nonlinear 
terms by additional changes of variables. 

For intermediate and strong noise, we obtained an approximation (I5.7P for the proba- 
bility ¥{N = 1} of spiking immediately, showing that the transition from rare to frequent 
spikes is governed by the distribution function $ of the normal law. Though we didn't 
obtain rigorous bounds on the principal eigenvalue and expected number of SAOs in this 
regime, simulations show a fairly good agreement with the approximation for P{A?" = 1}. 
The results on the Markov kernel contained in the appendix should in fact yield more 
precise information on Ao and the law of N, via approximations for the quasistationary 
distribution ttq. Generally speaking, however, we need better tools to approximate QSDs, 
principal eigenvalues and the spectral gap of substochastic Markov chains. 

Finally, let us note that the approach presented here should be applicable to other 
excitable systems involving oscillations. For instance, for some parameter values, the 
Morris-Lecar equations [ML81] admit a stable stationary point surrounded by an unstable 
and a stable periodic orbit. In a recent work [DGllj, Ditlevsen and Greenwood have 
combined this fact and results on linear oscillatory systems with noise [BG11] to relate 
the spike statistics to those of an integrate-and-fire model. It would be interesting to 
implement the Markov-chain approach in this situation as well. 

A Dynamics near the separatrix 

The appendix contains some of the more technical computations required for the proof 
of Theorem 14.21 We treat separately the dynamics near the separatrix, and during the 
remainder of an SAO. 

In this section, we use the equations in (£, z)-variables given by fj4.8f) to describe the 
dynamics in a neighbourhood of the separatrix. To be more specific, we will assume that 
z is small, of the order of some power of fj,, and that £ varies in an interval [-L, L], where 
the parameter L is given by ()4.15p . Let F± be the two broken lines defined in (14.131) 
and (|4.14p . Given an initial condition (—L,zq) G F^, our goal is to estimate where the 
sample path starting in (-L, Zq) hits F + for the first time. This will characterise the first 
part of the Markov kernel K. 
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A.l The linearised process 

Before analysing the full dynamics of (|4.8p we consider some approximations of the system. 
The fact that ^ ~ £ + t /2 for small z motivates the change of variable 



which transforms the system (|4.8p into 
du t = (-z t + 0(e)) dt + a x dW t {1) 



dz t = (& + tz t + 2u t z t + 0{ej) dt - a x t dW t {1) - 2a x u t dW t (l> + a 2 dW [ t 



(A.l) 



(A.2) 



where we write e = \/e(L 4 + cL 2 ). We choose an initial condition (0, Zq) at time to = — 2L. 
As a first approximation, consider the deterministic system 



du° = -z°dt, 

dz° = (/t + tz f °) dt . 



The solution of the second equation is given by 



y o _ J? 12 



z e 



-tjj/2 



to 



(A.3) 



(A.4) 



In particular, at time T = 2L, we have £t = L + ut — L and the location of the first-hitting 
point of F + is approximated by 



z T = Z() + fle TV2 f T e sV2 ds = Z0 + o{ ~l-j ) 
Jt 



(A.5) 



(A.6) 



As a second approximation, we incorporate the noise terms and consider the linear SDE 

du] = -z\ dt + a x dw\ X) , 
dz\ = (/} + tz}) dt - a x t dwl l) + a 2 dwi 2) . 

Let us now quantify the deviation between (u|,2^) and 
Proposition A.l. Let 



coo 



e *o + 



du 



to 



(A.7) 



T/ien there exists a constant M > suc/i t/iat /or all t ^ to, a// h,hi,h2 > and a// 

P e (o,/i 2 VM), 



sup J^Ljfi 



2(t - 1 ) r l/* 2 , __ 2 , 

exp<^ ---^(1 - Mp/i 7 y 



(A. 



and 



sup | tig — u® 



>hi + h2 I VC{s)ds 

Jtn 



^ 2exp 



2(t - to)fff 



to 

2(t-t ) 
H exp 



lh l(l-M P ^) 



8a 



(A.9) 
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Proof: The difference (x 1 ,^ 1 ) = (u 1 — u ,^ 1 — z°) satisfies the system 

dx] = -y\ dt + cri dW t {1) , 



dyj = t y] dt - a x t dW t {1) + a 2 dW i 



r(2) 



(A.10) 



The second equation admits the solution 

y \ = ~a 2 e * 2 /2 f e -s>/2 dW p _ ~ x e tV2 I 

J tn J t( 



se 



2 / 2 dWP =: 



1,1 . 1,2 

Vt +v t ■ 



(All) 



We first estimate . Let uq = to < u\ < ■ ■ ■ < uk = t be a partition of [to,t]. The 
Bernstein-like estimate [BG02[ Lemma 3.2] yields the bound 



t 



sup 



for any Hq > 0, where 

P k < exp 
The definition of £(s) implies 



a 2 



to 



K 



>H \ ^2^P k 



k=i 



2 CTg C(^fc) 



« 2 -s 2 



(A.12) 



(A.13) 



C(fl e «2_ s 2 = 1 _ _L f U " ^ u 2_ u 2 du ^ 1 _ e tg-„2 ^ _ (A. 14) 

C(«fe) COfc) Js Js 

Note that e*o = e 4 ^ 2 = 0(/i~ 27 ). For a uniform partition given by u k — Uk-i = P with 
p <C /x 27 , we can bound this last expression below by 



\-MpjX 



-2l 



for some constant M. This yields 
1,1 1 



sup 



[Us 



2(t - 1 ) 
P 



exp<j--^f (l-Mp/T 2 ^ 

Z (7 2 



(A.15) 



(A.16) 



1 2 

Doing the same for y s ' we obtain 



sup 



I 1,2 1 

\Vs I 



0^ss:t y / C(s) 



2(t - to) 



exp^-i§(l-Mp/i- 27 ) 

A (7 ^ 



(A.17) 



for any i2i > 0. Letting h = H$ + H\ with .ffo = H\ = h/2, we obtain (IA.8p . Now we can 
express x\ in terms of y\ by 



\ = - fy l s ds + d 1 fdWP 

Jto Jto 



Then the Bernstein inequality 



sup 



ax I dW. 

to 



(l) 



^ hi > ^ 2 exp 



2(t - t Q )a( 



yields (TO]) . 



(A.18) 

(A.19) 
□ 
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A. 2 The nonlinear equation 

We now turn to the analysis of the full system f|4.8[> . or, equivalently, (|A.2|) . Before that, 
we state a generalised Bernstein inequality that we will need several times in the sequel. 
Let Wt be an n-dimensional standard Brownian motion, and consider the martingale 

ft n rt 

M t = g(X s ,s) dW s = Y, 9i(X t ,t) dW t W , (A.20) 

J to i = l J to 

where g = (<7i,... ,g n ) takes values in M. n and the process X t is assumed to be adapted 
to the filtration generated by Wt- Then we have the following result (for the proof, 
see [BGK121 Lemma D.8]): 

Lemma A. 2. Assume that the integrand satisfies 

g(X t ,t)g(X t ,tf ^G(t) 2 (A.21) 
almost surely, for a deterministic function G(t), and that the integral 

V(t) = f G(s) 2 ds (A.22) 

Jtn 



is finite. Then 



for any x > 0. 



D ( sup M s > x) < e~ x2 / 2V W (A.23) 



Proposition A. 3. Assume zq = 0{fi l 7 ). There exist constants C,k,M > such that 
fort ^t^T + C(|log fi\~ l/2 ), all a < p, and H > 0, 



\zs~z° s \ ^ „\ ^CT ( f [H-M{T 2 fi 2 -^ + Te^)] i - 2/ - 2 
sup ' - sl ^ H } ^ — exp<^ -k± — - — = — > + e 

(A.24) 



and for all H' > 0, 



CT / f [#' - M(T 2 fi 2 ~^ + Tefi- 2 "/)] 



t ^t- j M 27 V I ^A- 27 



sup It*. - „2| > H> < ^ exp -. ^ r'"'* ;J + 



(A.25) 

Proof: The upper bound on t implies that e* 2//2 = 0(/2 -7 ). Thus it follows from (|A.3P 
and (fAT4|) that 



z° = ©(/i 1 - 7 ) and u° s = 0{Tfi 1 ^) (A.26) 
for to ^ s ^ i. Given h, h\,hi > 0, we introduce the stopping times 

r\ = infjs ^ i : \z\ - z° s \ ^ Ziy'C(s)} , 

r 2 = infja ^ to: \u\ - u° s \ >h x + h 2 J ^/C{s)ds^ . (A.27) 



The integral of \fCjs) is of order Tfi 7 at most. Thus choosing h = hi = hi = p, 
guarantees that 

z\ = 0{Ji l ' y ) and u\ = 0{Tfi l -~<) (A.28) 
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(D (A.30) 



for t ^ s ^ t A T\ A r 2 . For these values of h, hi and /i 2 , Proposition I A. II implies that 

F{n < t} < cTjl-^e-^ 2 /* 2 , 

IP{r 2 < t} < c27T 27 e -^ 2/ ^ 2 (A.29) 

for some constants k, c > 0. We consider the difference (x 2 ,y 2 ) = (ut, z t ) — {u\, z}), which 
satisfies the system of SDEs 

dx 2 = {-y 2 + 0{e))dt 

dy 2 = [ty 2 t + 2{u\ + x 2 )(zl + y 2 ) + 0(e)] dt - 2<r 1 (u t 1 + x 2 ) dW t (1) . 
We introduce a Lyapunov function Ut > defined by 

(U t - C ) 2 = (*' 2 ) 2 + fo 2 ) 2 . (A.31) 

The constant Co will be chosen in order to kill the second-order terms arising from Ito's 
formula. Let 

t* = inf{t ^ t : U t = 1} . (A.32) 
Applying Ito's formula and choosing Co of order ct 2 /! - ^ 1-7 ) yields 

dU t < [d + C a (*)l7 t ] dt + dWi 1 , (A.33) 

where (using the fact that a ^ p) 

C 1 = 0(T{i 2 - 2 ~<)+0{e) , 

C a (t) = t V0 + O(T/2 1-7 ) , (A.34) 
and <?(i) is at most of order 1 for t ^ t\ A t 2 A r*. Hence 

(•MtiAt 2 At* ptf\T\ AT2 At* 

^AnAraAr* < ^o+Cl(*-*o) + / C 2 (s)C/ s ds + c^ / </(«) dtf/ 1 . (A.35) 

J to J to 

We introduce a last stopping time 

r 3 = infjt > t : *i y <?0) d^ 1 > /i 3 | . (A.36) 
Then Lemma IA.2I implies 

F{r 3 < t} < e - * 3 ^/ 3 ? (A.37) 
for a K3 > 0. Applying Gronwall's lemma to (IA.35j) we get 

UtATiAT2Ar 3 AT* ^ [U t0 + Ci(t - t ) + /» 3 ] exp<^ / C 2 (u) du I 

= 0(T V~ 37 ) + O^T 2 ^ 1 - 27 ) + O^TpT 1 ) . (A.38) 

This shows in particular that r* > t, provided we take 7 small enough. Now (jAT24|) follows 
from the decomposition 

p( sup ^ i/l < p( sup ^ H — sup 



Jo<S<t \A>( S ) J Uo<s<tATiAT2AT3 yC(*j to^s^tAriAr 2 AT3 \/((s) . 

+ P{n > t} +P{r 2 > t} +P{r 3 > t} , (A.39) 
and (|A.25p is obtained in a similar way. □ 
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We can now derive bounds for the contribution of the motion near the separatrix to 
the Markov kernel. 

Proposition A. 4. Fix some 7 G (0, 1/4) and an initial condition (£0, zq) = (-L, zq) G F_ 
with \z \ = ©(/i 1 " 7 ). 

1. Assume c = 0. Then there exist constants C,K\,ho > such that the sample path 
starting in (£o>zo) will hit F+ for the first time at a point (£1,2:1) such that 

f[ Zi ^4 -/*} < ^expj-Ki^j . (A.40) 

2. If c ^ 0, but y/e ^ / u 1 + 2 7+^ f or som e 9 > 0, i/ien i/ie first-hitting point of F + always 
satisfies 

P{ki-4l ^ ^expj-Ki^} • (A.41) 

Proof: Consider first the case y/e ^ fi 1+2l+d . For any h > we can write 

P{|z T -4| ^P{|z T -4| ^/2-/i/2 2 - 27 }+P{|4-4| ^/i/i 2 - 27 } . (A.42) 

The first term on the right-hand side can be bounded, using (|A.24j) . by a term of order 
^-27 Q—Kifi /o" _ xhe conditions on 7 and yfe ensure that the error terms in the exponent 
in (|A.24j) are negligible. 

To bound the second term on the right-hand side, we note that (|A.3P implies that 
\z® — 4l ^ as or der /2 1_7 |r — T\. Furthermore, the definitions of r and u° imply that 
\ T -T\ = 2\u t - u° t \ + 0{jl 1 -^). This shows that 

P{ \z° - 4| > hfi 2 ~ 2 ^} < P{ |tit - u°| > /ici/i 1 " 7 - C2/U 1 " 7 } (A.43) 

for some constants 01,02 > 0. Taking /i = 2c2/c\ and using (|A.25|) yields a similar bound 
as for the first term. 

In the case c = 0, we can conclude in the same way by observing that zt is bounded 
below by its value for e = 0, the e-dependent term of dzt in (j4.8j) being positive. Thus we 
need no condition on \fe for the error terms in the exponent to be negligible. □ 

B Dynamics during an SAO 

B.l Action-angle-type variables 

In this section, we construct another set of coordinates allowing to describe the dynamics 
during a small-amplitude oscillation. Recall that in the limit e — > and jj, — >■ 0, the 
deterministic system (|4.4|) admits a first integral 

Q = 2ze- 2z - 2 ^ 2+1 . (B.l) 

The separatrix is given in this limit by Q = 0, while Q = 1 corresponds to the stationary 
point P. When e and ft are positive, we obtain 



Q 



2/2(1 - 2z) + Ve( ^|£ 4 + c[4z£ 2 - 6£ 2 + 8z 2 - 4z + l]^ 



e- 2 *-^ 1 . (B.2) 



Observe that if c = 0, the term of order yfe is strictly positive. 
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-2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 



Figure 11. Graph of u i-> /(u). 



In order to analyse the dynamics in more detail, it is useful to introduce an angle 
variable <fi. We define a coordinate transformation from (0, 1] x S 1 to R x R + by 



logQ 



■ sin ( 



3 =2 1+/ 



logQ 



COS « 



Here / : R — > (— 1, +oo) is defined as the solution of 

log(l + /(u))-/(u) = -2u 2 

such that 

sign /(it) = signu . 
The graph of / is plotted in Figure [TT1 

Lemma B.l. The function f has the following properties: 

• Lower bounds: 

f(u) > -1 and f{u) > 2u Vit G R . 

• Upper bounds: There exist constants C\,Ci > and a function r 
< r(u) < Cie" 1 " 2 " 2 , such that 



f(u) < C 2 u + 2u 2 
/(«) = 

Derivatives: f G C°° and 



l + e-^fl + r^)] 



Vu > 
Vu < 



/'(u) = 4u 



1 + 

/(«) 



/"(«) = 4 



1 + 

/(«) 



1-4 



u 



(B.3) 



(B.4) 
(B.5) 



(B.6) 
, with 

(B.7) 
(B.8) 



(B.9) 
(B.10) 
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There exists a constant M > such that 

< f"{u) s£ M VmG 



(B.ll) 



Proof: The results follow directly from the implicit function theorem and elementary 
calculus. □ 

We can now derive an expression for the SDE in coordinates (Q, (ft). To ease notation, 
we introduce the function 



X = X(Q, 
a parameter a > defined by 



logQ 



■ COS ( 



^ 2 ~ 2 I ~2 



(B.12) 
(B.13) 



and the two-dimensional Brownian motion dWt = (dW^ ,dW^ ) T . 

Proposition B.2. For z > 0, the system of SDEs fj4.8f) is equivalent to the system 

dQt = £/i(Qt, <ftt) dt + etyxCQt, t ) dW t 
d& = / 2 (Qt, dt + aMQt, <k) dW t , 



(B.14) 



where we introduced the following notations. 
• The new drift terms are of the form 



HQ, 

/ 2 (Q, 



-2Q 



a 1 



i + ^r q , £ {qa) + —RqAQi 
fj, /i 



fix) 

2X 



1 + 



2/itani 



log Q(l + /(*)) 



The remainders in the drift terms are bounded as follows. Let 

p{QA) 



VIlogQI if cos (ft ^ , 



(B.15) 
(B.16) 

(B.17) 



Then there exists a constant M\ > such that for all Q £ (0, 1) and all (ft E § 1 , 

^ MxllogQI 2 , \RqAQM ^ M lP {Q,cft) , 

\R^e(QA)\ < MillogQl^pCQ,^) , \R^{Q,<t>)\ ^ M lP (Q,(ft) 2 /\logQ\ . (B.18) 

Furthermore, if c = t/ien —f(X)RQ t£ (Q,(ft) ^ 0. 
TTie diffusion coefficients are given by 



tfti(Q,(ft) = [2V2^Q 



V-^gQ 



1 + f(X) 

2 l + /(X)cos0 or 2 1 /(X) 



tan< 



a v -logQ [l + /(X)]cos0' a logQl + /(X) 
There exists a constant M 2 > sitc/i i/iai /or a// Q £ (0, 1) and a// ^ £ S 1 , 

p(Q,0) 2 



t(Q, 



^M 2 Q z p(Q,, 



>XQ, 



logQP 



(B.19) 



(B.20) 
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Proof: The result follows from Ito's formula, by a straightforward though lengthy com- 
putation. The difference between the bounds obtained for cos (ft and cos (ft < is due 
to the fact that terms such as /(X)tan(0)/(1 + f(X)) can be bounded by a constant 
times \J — log Q in the first case, and by a constant times Q~ cos ^ in the second one, as a 
consequence of Lemma IB. 11 The fact that —f(X)RQ t£ is positive if c = follows from the 
positivity of the term of order y/e in (|B.2|) . □ 



B.2 Averaging 

In System (|B.14p . the variable Q changes more slowly than the variable (ft, which is a 
consequence of the fact that Q is a first integral when p, = e = a = 0. This suggests to 
use an averaging approach to analyse the dynamics. However, since the behaviour near 
(ft = 7r has already been considered in the previous section, using (£, z)-coordinates, we 
only need to consider (ft S [(fto, (fti], where — n < (fto < < (ft\ < n. 
We look for a change of variables of the form 

Q = Q + fiw(Q, (ft) (B.21) 

which eliminates the term of order fx in dQt- Ito's formula yields 

dw , , dw If d 2 w 2 d 2 w „ d 2 w 



dQ t = dQ t + fi— d(ft t + fi^ dQ t + -fi ( 7j7^ dQt + ^j^7 &Qt d<ftt + jr^ d(fti 



(B.22) 

dQt = Mfi + + 0(A) + 0(a 2 )) dt + af^+fi + ) dW t . (B.23) 



Replacing dQ t et deftt by their expressions in (IB.14j) . we get 

— nhr, -+- 

Thus choosing the function w in such a way that 

A + §£/a = (B.24) 

will decrease the order of the drift term in (|B.23p . We thus define the function w by the 
integral 

w(Q,0 = - l*¥§nr^i (B.25) 

which is well-defined (i.e., there are no resonances), since (|B.6j) shows that f2(Q><ft) is 
bounded below by a positive constant, for sufficiently small fl, e and a. 

Lemma B.3. Let (fto S (— n, — vr/2) and 0i € (7r/2,7r) 6e swc/i f/iai cos 2 (^o)j cos 2 (^i) b 
for some b G (0, 1). T/ien 



w(Q,(fti) X 



,2 a. ^ c ;„2. 



nsm' 0o n sin 01 

+ ^— -+ri(Q) 



sin (po sin <pi 



v^fogQ 
where the remainder terms satisfy 

r 1 (Q) = 0(Qlog(|lo g Q|)) , 



(l + r 2 (Q) + r e (Q)) , (B.26) 



Q- b Q~ 2b 



I log Q I V A |log Q| 

re(Q)^V~eO^^ + Q- b \lo g Q\^ , (B.27) 
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and r £ {Q) ^ is c = 0. Furthermore, the derivatives of w satisfy the bounds 

«? (e '* )=0 (7Kl)' ^«"> = ^'v^) ■ (B.28) 

and 

(B.29) 



sg 2 v^i^Qiy ' aw V^i^k 



Proof: We split the integral into three parts. Using the change of variables t = sin (ft 
and a partial fraction decomposition, we find that the leading part of the integral on 
[— 7r/2,7r/2] satisfies 

-/» 4QX(Q,#) „ / log(|logg|)N 



_„ /2 l+/(X(Q,<f.)) V 



Next we consider the integral on [fio, — 7r/2], The change of variables u = \J—2 log Q sin 0, 
(|B.8p and asymptotic properties of the error function imply 



-*' 2 aqx J± _ 2Q r v_2l0gQ du 



1 + /P0 



sin<W-21ogQ 1 + J^zM _ 1^ 



/•-V-21ogg 

2e / e~" /2 [l + Q(Qe-" / 2 )] dn (B.31) 

sin (poV— 21og<y 

O sin2 ^° r / 1 \ 2 i 

r°(M) +o(r > 



-2e 



v / -21ogQ(— sin< 



The integral on [tt/2, 0j] can be computed in a similar way. This yields the leading 
term in (|B.26|) . and the form of the remainders follows from (|B.18P with p = Q~ b . The 
bound on dw/d<p follows directly from ()B.24|) . while the bound on dw/dQ is obtained by 
computing the derivative of f\j fi- The bounds on second derivatives follow by similar 
computations. □ 

Notice that for the remainder f2(Q) to be small, we need that Q b 3> j2 and Q b S> cr 2 /A- 
Then the term r £ (Q) is of order y^ellog /i| 2 //i, which is small for /2/|log/2| 2 S> \fe. If 
that is the case, then w(Q,4>i) has order Q x ~ h ' / yjlog Q\. Otherwise, w(Q,4>i) has order 
\/eQ 1 ~ b \log Q\ 3 / 2 Ifi. In the sequel, we will sometimes bound l/yjlog Q\ by 1 to get simpler 
expressions. 



B.3 Computation of the kernel 

We can now proceed to the computation of the rotational part of the kernel of the Markov 
chain. Recall the broken lines F± introduced in (|4.13p and (|4.14p . For an initial condition 
(L,zq) £ F + , we want to compute the coordinates of the point (£, T ,z T ) at the first time 

T = w£{t>0:(£t,Zt)€F-} (B.32) 

that the path starting in (L,zq) hits F_. 
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We will assume that there is a /3 E (0, 1] such that 



{c-p)P ^ Z Q ^ z max < 



1 



(B.33) 



The (Q, </>)-coordinates of the initial condition are given by 

Qo = 2z e 1 " 2 ^ e" 2i2 > 2(c_ / u) /3+ ^ 



2L 2 



sin 



7 



- log Qo /3 + 7 

with cj)Q S (— 7r, —it/2). Thus Lemma IB. 31 applies with b = /3/(/3 + 7) < 1. Notice that 



(B.34) 



(B.35) 



Proposition B.4. Assume zq satisfies ()B,33p /or a /3 < 1. T/ien i/iere exists a constant 
k > suc/i i/iai i/ie following holds for sufficiently small \x and a . 

1. If \fe ^ /i/|log/2| 2 ; then with probability greater or equal than 



1 - e 



-nfi 2 /a 2 



zt) hits F_ for the first time at a point (— L, z\) such that 



zi = z + p,A(z ) + 



zq 



1-2*0 

TTie function A(zq) is given by 

p 2zo 



aV(z ) + O(~^ 1 -^+a 2 ^) 



A(z ) 



L(l - 2z ) 

and V(zq) is a random variable satisfying 
P{a\V(z )\ ^ h] s$ 2exp 



l + O(z log|logA|) +0 



^ 2 /2 2 ^ 



llog^ol 



V/i > 



(B.36) 



(B.37) 



(B.38) 



(B.39) 



2. If c = and -^/e > /i/|log/2| 2 ; t/ien (£t,2t) /ufcs i 7 - for the first time either at a point 
(—L,zi) such that Z\ is greater or equal than the right-hand side of (|B.37p . or at a 
point (£1, 1/2) with —L ^ £1 ^ 0, again with a probability bounded below by (|B.36j) . 

Proof: We first consider the case yfe ^ /2/|log/i| 2 . 

• Step 1 : To be able to bound various error terms, we need to assume that Qt stays 
bounded below. We thus introduce a second stopping time 



n = inf{t > 0: cos0 t < 0, Q™ 8 04 < (c_/i)^} 



(B.40) 



We start by showing that tAti is bounded with high probability. Proposition IB. 21 
implies the existence of a constant C > such that 



d& > Cdt + aMQt,<Pt)dW t 
Integrating this relation between and t, we get 



^0o + Ct + a [ MQs,0s)dW s 
Jo 



(B.41) 



(B.42) 
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Lemma IA.2I and (|B.20p provide the bound 

»tAn 



" i 4>2{Q s ,<Ps)<iw s 





^ h > ^ exp 



o- 



for some k > 0. Since by definition, 



j tAti 



2-K + h 
t A t\ > — — — ) ^ exp 



< 2ir, we get 



(B.43) 



(B.44) 



From now on, we work on the set f2i = {t A T\ ^ (2tt + 1)/C}, which has probability 

greater or equal 1 — e~ Kfl > a . 

Step 2 : The SDE (lB~22jl for Q 4 can be written 



dQt = QtKQt, <t>t) dt + aQMQt, <t>t) dW t , 
where the bounds in Proposition IB, 21 and Lemma IB .31 yield 

/(Qi) = 0(^)+aV- W ). 

MQ,4>)\\ 2 = o{^) . 

By Ito's formula, the variable Zt = log Qt satisfies 

dZ t = f(Z t , <p t ) dt + *$(Z t , <k) dW t , 
where f(Z, 4>) = f{e z , 0) + 0{a 2 lT 2 ^) and $(Z, <j>) = 4>{e z , <p). Setting 



(B.45) 



V t = / ^{Z s ,<j> 8 )dW 8 , 
Jo 

we obtain, integrating (|B.47|) and using the fact that /i 1 ~ 3/3 ^ A~ 2 ^, 

Z t = Z + dV + O{^ 1 -^+a 2 fl-^) . 
Another application of Lemma IA.2I yields 



a 2 



(B.46) 



(B.47) 



(B.48) 



(B.49) 



(B.50) 



for some k\ > 0. A convenient choice is h\ = jl 1 13 . From now on, we work on the set 
fii n n 2 , where ft 2 = {aV tATl < fl 1 '^} satisfies P(0 2 ) > 1 - e^ 1 ^ 2 . 
Step 3 : Returning to the variable Q, we get 



Qt = Qo e 



l + Otu^-^+a 2 ^ 



and thus 



Qt = Qoe aVt 



1 + 0(^2(1-/3) + ~2 --2^1 _^(Q tj 



(B.51) 



(B.52) 



Using the implicit function theorem and the upper bound on w, we get the a priori 
bound 

|9t ~ Qo1 = 0{j£~ p + a 2 /T^) . (B.53) 
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Step 4 : The a priori estimate ()B.53|) implies that on Oj n the sample path 
cannot hit F!_ on the part {— L ^ £ ^ 0, z = 1/2}. Indeed, this would imply that 
Q T ^> (c_/2) 7 , while Qq ^ a(c„/2) 7 with a = 2z max e 1_22max < 1. As a consequence, we 
would have (Q T — Qo)/Qo > (1 — a) /a, contradicting (|B.53|) . 

Let us now show that we also have t\ ^ r on Vt\ n ^2- Assume by contradiction that 
T\ < t. Then we have Q T1 = (c-.fl) 13 and cos 2 cj) T1 < /3/(/3 + 7), so that Q T1 = e>(/2^ +7 ). 
Thus jjLw(Q T1 ,4) T1 ) = 0(Q T1 ^ 1 - /3 ) = o{jl 1+r ). Together with the lower bound (IR341) 
on Qo, this implies that the right-hand side of ()B.52p is larger than a constant times 

at time t = t\. But this contradicts the fact that Q Tl = e>(//^ +7 ). 
Step 5 : The previous step implies that £ T = — L on £li n ^2- We can thus write 



4>t = g(Q T ) where sm(g(Q)) 
with g(Q) G (7r/2,7r). Notice that <?(Qo 



logQ L ' 

0. Furthermore, we have 



L 



V2Q(-logQ)3/2 C os( 5 (Q)) ' 
and thus Q t 9'{Qt) = C(l/|log /i| ) . Using this in the Taylor expansion 



(B.54) 



(B.55) 



w(Q T ,(j) T ) = W(Q , -</>o) + (Q T ~ Qo) 



^(Q0,g(Q0)) + ^(Qe,g(Qe)W(Q6 



which holds for some Qq £ (Qo,Q T ), yields the estimate 



W{Q T ,4> T ) w(Q ,-(po) 



Qo Qo 

Substitution in ()B.52p yields the more precise estimate 



Qi — Qc 



aV T ~w(Qo,-4>o) 



-A* 



Qc 



+ 0(^2(1-/3) + ~ 2-1- 
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(B.56) 
(B.57) 

(B.58) 



Step 6 : Finally, we return to the variable z\ = z T . Eliminating from the equa- 
tions (|R3|) . it can be expressed in terms of Q T as 



zi = G(Q T ) :-- 



1 



1 + / 



logQ 7 



L 2 



(B.59) 



Note that G(Qq) = Zq, while 

1 1 + f 

G'(Qo) 



log Qo _ J2 



zo 



2Qo 



f 



log Qo _ j2 



Qo(l - 2z ) 



(B.60) 



and 



G"{Q) 



1 + f 



logQ _ J2 



2Q 2 



logQ _ jjl 



logQ _ J2 



(B.61) 
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which has order Zq/Qq. The Taylor expansion 



zi = G(Qq) + (Qi - Qo)G'(Qo) + 



(Qi - Qo) 2 ^// 



G {Qo) 



(B.62) 



thus becomes 



. Qi — Qo zo . /0 
zi = zq-\ h O 



Qo 1 - 2z 



Qi — Q o 
Qo 



-zo 



By pT58|) . we have 

Qi — Qo 
Qo 

and Lemma lB.31 yields 

w(Qo,-<l>o) 2e 



Qo 



Qo 



L 



Q\ 



- cos <f>0 



+ 0(log|logQ c 



1 + 



I log Qc 



(B.63) 



(B.64) 



(B.65) 



Now (lBT34]l implies Q cos2 *° = e 22 °(2ezo) _1 and ci|logz | |logQ | ^ c 2 |log/2|. 
This completes the proof of the case y/e ^ /i/|log/2| 2 . 

In the case \fe > /i/|log fi\ 2 , we just use the fact that Qt is bounded below by its value in 
the previous case, as a consequence of ()B.28p . □ 

Corollary B.5. Assume that either c = or \fe ^ fi/\\og jl\ 2 . There exists a K2 > such 
that for an initial condition {L,z$) £ F- with zq ^ (c-/!) 1-7 , the first hitting of F+ occurs 

~ 2 I - 1 

at a height z\ ^ zq with probability larger than 1 — e~ K2 ^ /<T . 



Proof: It suffices to apply the previous result with j3 = 1 — 7 and h of order /2 7 . 



□ 
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